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ABSTRACT 


The  conventional  method  of  imposing  time  dependent  boundaury  condltioas  for  Kunge- 
Kutta  (RK)  time  advancement  reduces  the  formal  accuracy  of  the  space-time  method  to  first 
order  locally,  and  second  order  globally,  independently  of  the  spatial  operator.  This  counter 
intuitive  result  is  analyzed  in  this  paper. 

Two  methods  of  eliminating  this  problem  are  proposed  for  the  linear  constant  coefficient 
case:  1)  impose  the  exact  boundary  condition  only  at  the  end  of  the  complete  RK  cycle, 
2)  impose  consistent  intermediate  boundary  conditions  derived  from  the  physical  boundary 
condition  and  its  derivatives.  The  first  method,  while  retaining  the  RK  accuraicy  in  all 
cases,  results  in  a  scheme  with  much  reduced  CFL  condition,  rendering  the  RK  scheme 
less  attractive.  The  second  method  retains  the  same  allowable  time  step  as  the  periodic 
problem.  However  it  is  a  general  remedy  only  for  the  linear  case.  For  non-linear  hyperbolic 
equations  the  second  method  is  effective  only  for  for  RK  schemes  of  third  order  accuracy  or 
less.  Numerical  studies  are  presented  to  verify  the  efficacy  of  each  approach. 


’This  research  was  supported  by  the  National  Aeronautics  and  Space  Administration  under  NASA  Con¬ 
tract  No.  NASl-19480  while  the  authors  were  in  residence  at  the  Institute  for  Computer  Applications  in 
Science  and  Engineering  (ICASE),  NASA  Langley  Research  Center,  Hampton,  VA  23681. 


Introduction 

The  recent  interest  in  long  time  integration  is  due  to  the  need  to  tackle  problems  in  areas 
such  as  aero-acoustics,  electro-magnetics  and  others.  Tins  in  tmn  necessitates  working  with 
higher  order  accurate  spatial  differencing  operators.  In  many  cases  the  time-stepping  algo¬ 
rithm  of  choice  is  a  multi-stage  Runge-Kutta  (RK)  of  temporal  order  of  accuracy  comparable 
to  the  spatial  one. 

Several  bothersome  issues  arise  when  using  RK  methods  for  long  time  integration.  A 
principle  concern  is  the  imposition  of  numerical  botmdary  conditions  which  retain  the  formal 
accuracy  of  the  numerical  method  and  guarantee  numerical  stability  of  the  solution.  For 
example  (see  Trefethen  [1],  or  Carpenter  et.  al  [2],  [3])  Lax  and  GKS  stability  are  not 
sufficient  to  assure  that  the  there  is  no  exponentially  growing  temporal  error  when  using 
realistic  grids.  To  alleviate  this  “anomaly”  the  use  of  spatial  operators  which  have  specific 
semi-discrete  energy  norms  has  been  proposed  [4],  [5],  [6],  [7].  These  papers  have  primarily 
focused  on  the  semi-discrete  form  of  the  equations. 

The  present  paper,  however,  deals  with  the  loss  of  accuracy  due  to  the  imposition  of 
time  dependent  bovmdary  conditions  g{t)^  dictated  by  the  physics  of  the  problem.  The 
conventional  way  of  dealing  with  the  uncertainty  of  what  is  happening  at  the  intermediate 
stages  of  the  RK  time  advancement  is  to  impose  at  t  -4-  the  boundary  value  g{t  -f  ot„5t), 
where  is  the  coefficient  appropriate  to  the  particular  stage  of  the  given  RK  algorithm. 

It  will  be  shown  that  this  conventional  boundary  condition  imposition  leads  to  a  numerical 
scheme  which  is  only  first  order  accurate  in  the  naghborhood  of  the  boundary,  leading  to 
a  global  accuracy  of  second  order  only.  Another  approach  is  to  treat  the  time-dependent 
boundary  .condition,  g{t),  as  a  source  term  in  the  governing  partial  differential  equation 
(p.d.e),  thereby  avoiding  the  need  to  formally  specify  intermediate  boimdary  conditions. 
However,  it  can  be  shown  that  procedure  is  equivalent  to  the  conventional  method  with  its 
attended  problems  (see  Appendix  A).  A  third  natural  procedure  is  indeed  not  to  specify  any 
intermediate  boundary  condition.  ()ut^ ,  |o  obtain  them  from  the  inner  scheme.  This  method 
retains  the  accuracy  of  the  spatial  operator,  but  significantly  reduces  the  allowable  time  step 
for  stability,  rendering  the  scheme  less  attractive. 

In  section  2)  we  analyze  and  pinpoint  the  reasons  for  the  deterioration  of  the  accuracy  and 
provide  a  simple  recipe’  for  restoring  the  accuracy  in  the  case  of  linear,  constant  coefficient 
hyperbolic  system  of  p.d.e.’s. 

In  section  3)  we  deal  with  a  non-linear  hyperbolic  system  of  conservation  laws.  The  reme¬ 
dies  provided  in  Sections  2)  and  3)  can  not  be  generalized  to  non-linear  systems  integrated 
by  RK  schemes  of  arbitrary  order  of  accuracy.  We  show  that  for  the  RK  of  third  order,  the 
remedy  of  section  2)  is  effective. 
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2*  The  linear  caie 


To  illustrate  the  phenomenon  of  loss  of  accuracy  due  to  the  conventional  imposition  of 
inflow  boundary  conditions,  we  consider  the  following  problem: 

(1) 
(2) 


dv  ■ 

~  -  (QV’(i)).  I  =  i  >  0  (3) 

=  5(0t  (‘"I) 

where  V  =  Vi^  ~  [uo  is  the  semi-discrete  approximation  w'hich  converges  to  u{xi^i) 

at  the  spatial  grid  points  X;  (for  stable  discretizations);  and  Q  is  the  differentiation  matrix 
representation  of  the  derivative  operator  The  specific  form  of  Q  depends  on  the  algorithm 
used  and  in  particular  on  the  order  of  accuracy.  For  all  finite  difference  operators  on  uniform 
grids  (which  suffices  for  the  present  purpose  of  illustration)  we  may  write  C}  ”  where 
h  is  the  mesh  spacing. 

The  demonstration  of  accuracy  deterioration  will  be  shown  for  the  four-stage  '‘classical” 
RK  scheme,  which  is  one  of  the  most  common  RK  time  advancing  schemes.  For  the  analysis 
to  make  sense  we  assume  that  the  differentiation  matrix  is  at  least  of  fourth-order  accuracy 
up  to  the  boundary.  It  should  be  noted,  however,  that  thi.‘5  illustration  could  be  carried  out 
for  any  RK  algorithm. 

The  above  mentioned  four-stage  integration,  together  with  the  conventionally  imposed 
boundary  conditions,  is  implemented  as  follows: 

+  5(DV'”). 

^0  =  9{i  +  -2  ) 

u?  =  ur  + 

V 

uf  -  u."  -t-  X{DV% 
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1  <  I  <  (5) 

1  <  i  <  A'  (6) 

1  <  t  <  .V  (7) 


4  =  9{t  +  it) 

^  vP  ^  jl(DV^)i^2(DV^)i  +  2(DV^)i  +  (DV^)i]  1<  i  <N  (8) 

=  9(t  +  it) 

where  A  =  ^. 

Equations  (5)  -  (8)  take  the  semi-discrete  variable  from  the  time  level  n,  to  Vi(ii-it) 
at  time  level  n  4-  f . 

For  the  purpose  of  analysis,  the  above  system  is  rewritten  in  the  following  form,  again 
with  V  ~  [uo, 


l/>  =  4-  ^£)V»  -f.  G^eo  (S) 

1/2  =  V"  +  eo  (10) 

+  ADl/2  4-  6?^  Co  (11) 

V/^+i  =  \/“  +  4-  2DV^  +  2DV^  4-  eo  (12) 

0 

where  eo  =  [1,0,  ...0]  ,  and 

ff>  =  <;{(  +  |)  -  <  -  icDV").  (13) 

O'  -  J(i  +  |)  -  vS  -  ^(DV'h  (14) 

O'  =  S((  +  6i)  -  <  -  ACOr^jo  (15) 

G=  =  9((  i-  St)  -  uj  -  -lOV"  +  2DV‘  +  2DV'  +  CV''’]o  (16) 
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Substitution  of  (9)  into  (10)  and  the  result,  into  (11),  etc.  leads  to  the  following  expression 
for  1/"+': 


K-‘+'  =  [/  4-  4d  4-  t 

l!  2!  3!  4! 

+  G«(jC  +  jD'  + 

+  G'l^D  +  jD'y 

O.  /::2fAj_)i^o 

■  -  tQ-  t 

4-  G^e°  (17) 
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To  determine  the  formal  accuracy  of  (17)  we  substitute  the  exact  solution  (pro¬ 

jected  at  the  points  cCj),  for  i;-'.  Under  the  assumption  on  the  order  of  the  differentiation 
matrix  Qj  it  is  clear  that  D^V'^  in  equation  (17)  can  be  replaced,  up  to  fourth-order  accuracy, 
by  ^*[(“1)^^'^!^]  at  the  points  Equation  (17)  then  becomes 

=  u{i  +  6t)  +  Oi[Stf,  [6xf) 

+  G^i^D  +  jD^  +  ~D^]e° 

-f  G^[^D]e° 

-r  (IH) 

Applying  the  same  procedures  to  equations  (13)  -  (16)  we  get; 


<3°  =  +  |)  -  9(1)  -  §9'(l) 

G"  =  9(<  +  f )  -  9(0  -  ^9'(<)  -  ^9"(0  -  G“ 

=  g(t  +  ft)  -  j(0  -  (ft)9'(0  -  ^9"(0  -  ^9"'(0  -  A4  G>  -  G” 

G’  =  f,((  +  ft)  -  s{i)  ~  (St)g\t)  -  -  i^9"'(0  “  ^9""(0 

“  IgdoO  +  “^^00  +  -Y^«OoJ^ 

"  [3^00  +  ^^^00]  G' 

~  ^^00  G^  (19) 

where  dg^  =  (Deo)c,  dgg  —  (D^eo)o,  dig  =  (Z)^eo)o,  and  (77^60)0  is  the  first  element  of 
the  vector  D^eo\  1  <  <  3.  A  Taylor  series  expansion  of  equation  (19)  clearly  shows  that 

G°  =  (9([(5t]^);  as  it  does  for  G',  G^  and  G^,  for  arbitrary  A,  In  addition,  it  can  be  shown 
that  the  vectors  D'^eo  (1  <  7:  <  3)  are  linearly  independent.  Substituting  equation  (19)  into 
equation  (18)  yields  the  expression 


{D^eo)g"  +  O{{6tf^) 


(■20; 
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Rearranging  equation  (20)  we  have 


yn+l 


(2i) 


which  means  that  the  RK  approximation  is  locally  (near  the  boundary)  only  first  order 
accurate.  The  locality  is  assured  by  the  finite  support  of  D  for  explicit  schemes,  and  by  geo¬ 
metrically  decaying  coefficients  in  space  for  compact  schemes.  Hov/ever,  the  overall  accuracy 
in  tlie  Linf  norm  cannot  exceed  second  order  (see  Gustafsson  [8]). 

The  remedy  for  the  dilemma  posed  by  equation  (21)  suggests  itself  when  one  examines 
equa.tions  (17)  and  (18).  We  note  that  if  we  set  —  G^  =■  0,  then  V'"'"*'’  = 

u(t  4-dO -t- 0([<5^j*)  -h  e.oG^.  But,  with  =  G*  -  =  0.  then  G^  =  0([dt]®),  and  we 

have  tli('  rorrect  order  fur  V’’  ’*"'.  To  achieve  these  identities  v  e  specifically  use  in  equations 
(5)-(8)  the  following  expressions  for  the  intermediate  boundary  conditions; 


=  s{‘)  ->■ 

(22) 

(23) 

=  s(l)  +  (S‘)s'(‘}  +  ^sr"'(0 

(24) 

Note  that  equations  (22  -  24)  are  precisely  the  intermediate  boundary  conditions  obtained  by 
numerically  integrating  ^  w’ith  the  clcissical  fourth  order  RK  scheme.  Specifically, 

replacing  the  third  order  equation  in  u  with  the  system  of  first  order  equations  for  the 
unknown  functions  Ko,  vto  and  vtto. 


d(V)o 

(it 

=  (uf)o 

(25) 

d(vi)o 

(it 

(26) 

d{vtt)o 

dt 

=  9"'{i) 

(27) 

on  the  boundary,  and  integrating  with  the  standard  fourth  order  RK  scheme,  yields  preci.sely 
the  intermediate  boundary  conditions  proposed  in  equation  (24).  .Solving  the  third  order 
o.d.e  (27)  on  the  boundary,  is  a  remedy  that  can  be  applied  to  any  four  stage  fourth  order 
RK  scheme,  and  thus  provides  a  sinqde  and  general  means  of  implementing  the  correct  inter¬ 
mediate  l)oundary  conditions.  For  three  stage  third  order  RK  schemes,  solving 


o 


on  the  boundary  pro%'ides  a  genera!  technique  for  imposing  the  correct  boundary  conditions 
on  the  linear  problem. 

At  this  point  we  comment  on  why  the  above  predicted  phenomena  has  not  been  observed 
previously  by  the  practitioners  of  higher  order  methods.  From  equation  (21)  we  see  that 
the  leading  error  coefBcient  in  the  expression  is  proportional  to  (A)‘‘5a;.  For  example,  if  D 
represents  a  fourth  order  central  difference  operator  with  fourth  order  boundary  closures, 
and  u  =  then  the  error  at  the  point  next  to  the  boundary  becomes 

^  t(184320  -  1.520A)A'‘oV  ,  (-{•  1 72800  +  23040A_- 190A^J^^^ 

~  79626240  79626240 

Thus  if  A  <<  1  for  a  giv^n  Si  it  is  possible  that.  \‘*Sx  <  (<*'0”*'  However,  if  one  refines  t.ie 
grid,  or  ''onverscly,  if  (Uie  runs  llic  ("(Jiiiputalion  at  the  allowable  stalrie  CFL  (A  —  0(1)), 
then  the  first  order  error  becomes  apparent.  A  detailed  study  of  this  effect  is  now  presented. 

We  begin  by  showing  the  the  results  from  a  grid  refinement  study  using  three  high-order 
central  difference  schemes;  1)  (3, 3-4-3, 3):  a  fourth  order  explicit  interior  with  two  third-order 
explicit  boundary  points  at  each  end  of  the  domain,  2)  (5^,5^  ~  6  —  5^,5^)[2];  a  sixth-order 
compact  interior,  with  two  fifth-order  boundary  closures  at  each  end  of  the  domain,  and  3) 
(5,5,5,5,5,5-6-5,o, 5,5,5, 5)[6];  a  sixth-order  explicit,  with  six  fifth-order  boundary  points  at 
each  end  of  the  domain.  (See  Carpenter  et  al[2]  for  detail  of  the  high-order  nomenclature). 
Ail  schemes  are  GKS  and  time  stable  for  the  scalar  hyperbolic  equation.  In  addition,  when 
used  in  conjunction  with  the  SAT  [3j  boundary  procedure,  the  semi-discrete  form  of  scheme 
3)  can  be  shown  to  be  time  stable  for  the  constant  coefficient  hyperbolic  system.  An  explicit, 
a  compact  implicit,  and  a  scheme  which  satisfies  the  summation  by  parts  energy  norm  were 
chosen  to  illustrate  the  generality  of  the  problem,  eis  well  as  the  remedy. 

The  model  problem  used  to  test  the  schemes  is  the  scalar  hyperbolic  equation 


du  du 
dt  dx 

=  0,  0<x<l,  t>0 

(29) 

u(0,t) 

=  sin  27r(-“t),  t  >  0 

(30) 

u(x,0)  = 

sin 27r(x),  0  <  x  <  1 , 

(31) 

The  exact  solution  is 
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u{x^i)  «  siti  27r(s -- i),  0<ce<1,  t>0  (32) 

In  all  cases  the  solution  was  advanced  in  time  using  the  classical  fourth-order  RK  scheme 
with  boundary  conditions  imposed  as  in  equations  (5)  -  (8). 

Upon  grid  refinement  with  a  vanishingly  small  CFL,  all  schemes  recover  their  theoretical 
accuracy  of  fourth-order,  sixth-order  and  sixth-order,  respectively.  Table  (La)  shows  that 
all  schemes  converge  at  a  fourth-order  rate  for  moderate  resolutions,  but  gradually  degrade 
to  an  asymptotic  rate  of  2.5.  This  trend  is  representative  of  high-order  explicit  or  compact 
spatial  operators  advanced  in  time  with  a  fourth-order  RK  scheme,  with  or  without  the  SAT 
boundary  procedure. 


(3,3-4-3,3) 

— 

_ 

(5L  5^  -  6  —  5'^,  5^) 

(5,5,5.5..5-6-) 

CFL  =  2.0 

Conv 

CFL  =  1.5 

Conv 

Grid 

log  Li 

Rate 

log  Li 

Rate 

41 

-2.371 

-1.537 

81 

-3.570 

3.98 

-4.242 

4.02 

-2.221 

161 

-4.678 

3.68 

-5.450 

4.02 

-2.956 

321 

-5.593 

3.04 

-6.635 

3.94 

-3.704 

641 

-6.380 

2.61 

-7.711 

3.57 

-4.455 

2.49 

1281 

-7.138 

2.52 

-8.586 

2.91 

-5.208 

2.50 

Table  La:  Conventional  imposition  of  boundary  condition  n(0,t)  =  g{t),  for  three 

high-order  finite  difference  schemes. 


Two  possible  remedies  have  been  suggested  to  rectify  the  loss  of  accuracy  for  the  linear 
problem:  1)  do  not  impose  intermediate  physic  umdary  conditions,  and  2)  impose  con¬ 
sistent  intermediate  boundary  conditions  derive,  =  ;m  the  physical  boundary  conditions  and 
their  derivatives.  (Or  alternatively,  solve  the  related  system  of  equations  on  the  boundary). 

Not  imposing  intermediate  physical  boundary  conditions  results  in  a  fully  discrete  scheme 
which  is  formally  fourth  order  accurate.  (By  definition  a  fourth  order  scheme  in  space  and 
time  will  remain  fourth  order  in  the  absence  of  boundary  conditions).  A  problem  v/ith  this 
remedy,  however,  is  that  the  stability  of  the  scheme  is  greatly  reduced.  When  using  the  RK4 
scheme  with  various  finite  difference  operators,  at  legist  of  a  factor  of  two  (and  often  much 
more)  decrease  in  CFL  was  observed. 

The  alternative  remedy  is  to  impose  consistent  intermediate  boundary  conditions,  derived 
from  the  physical  boundary  conditions  and  its  derivatives.  For  the  scalar  advection  defined 
by  equations  (29)  -  (30)  it  is  sufficient  to  solve  the  derivative  boundary  conditions  described 
by  the  system  of  equations  (27). 


7 


Table  (l.b)  shows  the  results  of  a  grid  refinement  study  using  the  derivative  form  of  the 
boundary  condition. 


■a:aK.Lfi«iTOIiai 

Conv 

CFL  =  1.4 

Conv 

CFL  ^  1.5 

Conv 

Grid 

log  Li 

Rate 

log  Li 

Rate 

log  Li 

Rate 

-2.394^ 

-2.490 

81 

-3.613 

4.05 

-4.282 

3.97 

-3.767 

4.24 

161 

-4.817 

4.00 

-5.486 

3.99 

-5.076 

4.35 

321 

-6.019 

3.99 

-6.687 

3.99 

4.32 

641 

-7.222 

3.99 

-7.891 

4.00 

-7.655 

4.25 

1281 

-8.425 

4.00 

-9.099 

4.01 

-8.911 

4.17 

Table  l.b:  New  physical  boundary  condition  imposed  as  three 

high-order  finite  difference  schemes. 


Fourth-order  temporal  accuracy  is  recovered  for  all  methods  with  this  approach.  In  addition, 
the  same  maximum  CFL  was  achieved  in  all  cases,  as  was  possible  with  the  conventional 
boundary  conditions.  Similar  results  have  been  obtained  for  three  stage  third  order  and  six 
stage  fifth  order  RK  schemes  as  well. 

Section  3.  Non-Linear  Hyperbolic  System 

We  have  shown  that  the  traditional  imposition  of  time-dependent  boundary  conditions 
causes  a  degradation  of  the  formal  accuracy  to  first-order.  We  have  also  shov/n  how  to  elim¬ 
inate  the  problem,  given  the  exact  boundary  condition  and  its  derivatives  on  the  boundary. 
We  shall  now  show  that  boundary  treatment  similar  to  that  resulting  from  the  linear  analysis 
of  section  2,  is  also  valid  for  third  order  RK  schemes,  even  for  the  Ccise  of  a  system  of  non¬ 
linear  conservation  laws.  The  procedure  does  not  obviously  generalize  to  higher  order  RK 
integration,  however.  (The  only  technique  available  that  does  directly  extend  to  non-linear 
hyperbolic  systems  is  not  imposing  intermediate  physical  boundary  conditions.) 

Consider  the  system 


dt 

0{o,t) 


=  A{U)  ~  0  <  X  <  1;  f  >  0 

ax  ox 

=  m 


(33) 

(34) 


where  A{U)  is  the  .Jacobian  of  f{U)  with  respect  to  V .  We  perform  the  integration  using  the 
third  order  RK  scheme  attributed  to  Heun.  The  boundary  conditions  at  the  intermediate 
time  levels  arc  obtained  from  .solving  the  boundary  equation  =  ^'[t]  with  Heuii’s 

method. 
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vl 

-  w7 

+  j(o/(r")l.- 

VI 

<  N 

(35) 

= 

Si  ^ 

+  J9«, 

= 

+  jlOfiS')]: 

1  < 

i  <  N 

(36) 

Wo 

=  ff(i} 

+  -TrS  (i)  d* 

2{St)^ 

9 

•r(^) 

^+1 

=  W  +  jj;  + 

1  <  i  <N 

(37) 

-Vr+l 

Wq 

-  9{i 

4-  Si) 

Follov  ing  the  notation  of  equations  (9)  -  (12)  we  have- 


+  ^[Df{r)]i  +  &eo  (38) 

+  “(n/V)].'  +  eo  (39) 

-  iu7  4-  h-  +  eo  (40) 

'I  ^  4 

where,  it  is  clear  from  equation  (35)  -  (37)  that  here: 

G°  =  3(i)  +  f  s'!')  -  «o  -  (■“) 

5'  =  m  i-  ~m  +  ^9"(<)  -  j|e/(5>)|o  («) 

=  S(l^il)  -  jv^-  -  ji?J  -  jAiO/V)i„  (43) 


As  in  the  linear  case  described  by  equation  (18),  will  be  given  by  a  linear  combination 
of  (7°  and  O',  plus  terms  proportional  to  ((50“*-  We  now  show  from  equ..tions  (41)  and  (42) 
that  t?’  =  0,  and  G’  =  G([^t]^),  and  thus  the  overall  accuracy  cf  is  0([di]^). 

As  in  section  2,  in  order  to  obtain  the  truncation  error  we  substitute  u{x,J)  for  7-  The 
vector  G'°  immediately  follows  as 

=  9(1)  +  |,f (0  -  (ff(o  +  |/i(5) )»  m) 

=  9(1)  +  J^(‘)  -  (“(1)+ 
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ajid  io  a  similar  vain, 


(?'  =  m  +  ^5^(0  +  ^j"«)  -  { ff(«)  -  + 1 /.(iJ)]}"  (46) 

=  m  +  x9"(‘)  +  ^r(i)  •*■ 

(47) 

but  f^ui  -  A{u)ut  -  ft  and  ■§^[Ji\  ~  §i{L)  ~  Thus 


2(6ty‘ 

=  m  +  •y^(o  +  -v- 
G'  =  om?) 


rm  -  ( -  ^s,  + 


2lSt)\^ 


u»}o  d*  G{[6t]  ) 


(48) 


The  consequence  of  equations  (45)  and  (48),  substituted  into  equation  (40)  with  — ♦  u(a;,  t) 

is  that  the  error  is  proportional  to 


-  Ui(t  +  <?t)| 

St 


0([dtp);  t  =  0,l,...,m 


(49) 


near  the  boundary,  (m  finite),  and  proportional  to  (9([^t]'^)  in  the  interior.  In  other  words 
,  the  boundary  conditions  prescribed  in  equations  (35)  -  (37)  give  us  the  required  overall 
accuracy.  Unfortunately,  this  procedure  does  not  seem  to  generalize  to  RK  integrations  of 
order  higher  than  three  for  the  case  of  non-linear  systems.  The  main  reason  is  that  beyond 
U(,,  we  will  need  to  use  the  “Jacobian  of  the  Jacobian”  which  can  not  be  related  simply  to 
the  temporal  derivatives  of  the  vector  u{x,t). 


4.  Conclusion 

We  have  shown  that  the  imposition  of  inflow  boundary  conditions  at  the  intermediate 
steps  of  Runge-Kutta  algorithms  for  hyperbolic  P.D.E.’s  has  to  be  done  in  a  counter-intuitive 
manner,  if  oue  is  to  preserve  the  overall  accuracy  of  the  scheme.  The  conventional  (or 
"natural”)  way  of  assigning  at  level  n  +  a,,  (0  <  a,  <  1)  the  value  g[t  a^St)  degrades 
the  scheme  to  be  of  first  order  accuracy  near  the  boundary  and  of  second  order  accuracy 
overall.  We  presented  ways  to  deal  with  the  linear  case  of  general  order  and  with  systems  of 
conservation  laws  in  the  ca.se  of  third  order  RK  integration.  Much  work  remains  to  be  done 
for  the  non-linear  case  and  linear  problems  with  variable  coefficients. 
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Appendix  A 

Once  the  spatial  operator  has  been  chosen,  a  P.D.E.  becomes  a  system  of  O.D.E.’s,  plus 
a  boundary  term.  If  the  boundary  term  i?  treated  as  a  source  term,  then  entire  system  can 
then  be  treated  by  conventional  techniques.  We  show  that  this  technique  suffers  from  the 
same  loss  of  accuracy  at  the  boundary  as  was  discussed  earlier. 

We  start  with  the  governing  equations 


du  du 
dt  dx 
u(0,0 


0  0  <  a:  <  1,  f  >  0 

dii) 


I'hc  somi-di.scrctizcd  version  of  (AI.  1)  -  (AI.  2)  is 


(AI.  1) 
(AI.  2) 


^  t  =  (A1.3) 

Uo(f)  =  g{i)  (AI.  4) 


where  V  =  =  0,...,  A^  Define  the  reduced  identity  nratrix  P  by 


'  f  ''T' 

where  the  matrix  P  is  dimension  (N-M,N  +  1).  Similarly,  Cq  —  [1,0,  ...,0]  ,  and  has  dimension 
(N  +  1).  Next  define  the  new  variables  V  -  PV,  D  —  PDP  and  qo  =  PDe^.  Now  replace 
equation  (AI.  1)  by  the  fully  discrete  source  equations: 


= 

+  jKDV)? 

+  9{i)(go')j]  1  <  i 

VI 

(AI. 

5) 

= 

+  jKDV); 

+  s{^  +  y)(9o);] 

(AI, 

6) 

- 

+  .'((ov;; 

+  P(^  +  y)(7o)j] 

(AI, 

7) 

On+l 

; 

D 

,  +  2{DV'),  f  2(DK'), 

+  (DV^):  1 

+  'I'ldit]  +  4  5(1  +  y)  +  g{t  +  ^t)](qo)j 

(AI. 

8) 
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Rewriting  equation  (AI.  6)  in  terms  of  V  and  D  one  obtains: 


(PV')i  =  {PV')i  +  -liPDPV")!  +  g{t)PD(e,)i]  \  <j  <  N 

=  {PV'h  +  ^{PDV'h  =  P\{V")  +  (AI.  9) 

It  is  clear  that  the  reduced  vector  PV'  is  identical  to  the  vector  at  time  level  1  obtained 
from  the  conventional  imposition  of  boundary  conditions  for  1  <  j  <  N  [see  equation  (5)]. 
Noting  this  equivalence,  equation  (AI.  6)  can  now  be  interpreted  as 


P{V%  =  r^MPDPV'),  +  i,(i  +  |)PC(e„),;  1  <  j  <  N 

=  P(V'),  +  ~{PDV')i  =  Picn  +  j(£IV')Ll  (AI.  10) 

The  reduced  vector  PV^  is  identical  to  that  obtained  with  the  conventional  boundary  condi¬ 
tions  for  1  <  j  <  N,  This  procedure  can  be  used  to  show'  that  each  stage  of  the  conventional 
boundary  procedure  is  equivalent  to  that  obtained  from  solving  equations  (AI.  5)  -  (AI.  8). 
Thus,  the  entire  procedures  are  equivalent. 
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